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We solve the coherent multiple Andreev reflection (MAR) problem and calculate current-voltage characteris- 
tics (rVCs) for Josephson SINIS junctions, where S are local-equilibrium superconducting reservoirs, I denotes 
tunnel barriers, and N is a short diffusive normal wire, the length of which is much smaller than the coherence 
length, and the resistance is much smaller than the resistance of the tunnel barriers. The charge transport regime 
in such junctions qualitatively depends on a characteristic value 7 = tjA of relative phase shifts between the 
electrons and retro-reflected holes accumulated during the dwell time Tj. In the limit of small electron-hole de- 
phasing y <C 1, our solution recovers a known formula for a short mesoscopic connector extended to the MAR 
regime. At large dephasing, the subharmonic gap structure in the JVC scales with y^' , which thus plays the role 
of an effective tunneling parameter. In this limit, the even gap subharmonics are resonantly enhanced, and the 
IVC exhibits portions with negative differential resistance. 

PACS numbers: 74.50.+r, 74.45.+C 
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I. INTRODUCTION 

Accurate theoretical description of nonequilibrium charge 
transport in Josephson junctions is an important and active 
research field. The concept of Multiple Andreev Reflections 
(MAR)' is a universal framework explaining the nature of dis- 
sipative current in different types of junctions. The quasipar- 
ticles injected in the junction at applied voltage smaller than 
the superconducting energy gap, eV < 2A, can only escape 
into reservoirs at zero temperature by multiple traversing the 
junction due to repeated Andreev reflections, each time gain- 
ing energy eV. Such a process generates coherent transfer of 
multiple electron charge ne — 2A/V across the junction. Such 
a mechanism is important at small temperatures when the sin- 
gle particle tunneling is exponentially weak, and the Andreev 
transport becomes dominant. The most complete quantitative 
MAR theory has been developed for ballistic contacts. The 
central element here is the solution for a single conducting 
channel with given transmissivity, which enters the sum over 
conducting channels in the net cuiTent. Various approaches 
for constructing such a solution have been developed based on 
the scattering theory, tunneling Hamiltonian, and quasiclassi- 
cal Green's functions. ^"^ 

Extension of the theory has been suggested for disordered 
weak links, e.g., diffusive constrictions and superconductor- 
insulator-superconductor (SIS) tunnel junctions with disor- 
dered insulating layers. In this case the summation over the 
channels is performed by taking into account the distribution 
of random transmission eigenvalues for corresponding struc- 
ture in the normal state. ^"'^ This method, however, is only valid 
for short superconducting weak links, in which the dwell time, 
i.e., the time of quasiparticle diffusion along the whole MAR 
staircase, is small on a quantum time scale defined by the 
inverse order parameter A^'. In junctions with large dwell 
time, e.g., long superconductor-normal metal-superconductor 
(SNS) junctions and double-barrier SINIS junctions contain- 
ing Breit-Wigner resonances, the method does not work be- 
cause of strong energy dispersion of the scattering ampli- 
tudes. The electron- and hole-like quasiparticles acquire dif- 
ferent scattering phases during propagation through the junc- 



tion (electron-hole dephasing), which results in the single- 
channel MAR current being dependent on the scattering am- 
plitudes rather than the scattering coefficients.'''"'^ Thus the 
distribution of transmission eigenvalues becomes irrelevant 
and is to be replaced by the statistics of scattering amplitudes, 
which is generally unknown except of some particular cases. '— 

In order to investigate the MAR problem in diffusive junc- 
tions with large electron-hole dephasing, one has to directly 
solve quasiclassical Keldysh-Green's function equations.'— 
This task, however, is technically demanding because a non- 
linearity of equations and nonstationary character of the MAR 
problem lead to hardly tractable two-dimensional infinite 
chains of Green's function harmonics. So far the problem has 
been analytically solved for superconducting constrictions,'— 
and for the opposite case of very long SNS and SINIS junc- 
tions, where the Josephson cuiTent is completely suppressed 
and a incoherent MAR theory can be formulated'^ '^ analo- 
gous to the one for ballistic junctions. '** A perturb ative scheme 
for a superconducting film interrupted by a tunnel junction has 
been suggested,'^ however no general methods to analytically 
treat the problem in SNS junctions of intermediate lengths ex- 
ists up to now, and even numerical solution presents a difficult 
and time consuming task.^** 

In this paper we present a study of the junction for which 
the coherent MAR problem can be relatively easily solved and 
fully investigated both analytically and numerically for the en- 
tire range from weak to strong electron-hole dephasing. We 
consider a SINIS junction with opaque NIS interfaces having 
equal resistances R much larger than the resistance Rn of the 
diffusive normal wire, the length 2d of which is supposed to 
be much smaller than the coherence length -^/D/A (D is the 
diffusion coefficient, Ti= l).ln such junctions, the dwell time 

is characterized by the paramete r '^'^ '" — 



Y = T^/A = , 

Rn Eih 



(1) 



where Ejh = D/(2c/)^ is the Thouless energy. This parameter 
defines the amplitude of the Josephson current and the magni- 
tude of the minigap in the energy spectrum within the normal 
wire A^, ^ A/(l + 7). In short junctions, the Thouless energy 



2 



is large compared to A, therefore, the parameter 7 may have 
arbitrary value depending on the junction length and trans- 
parency. Correspondingly, the value of the minigap may vary 
between and A reflecting the crossover from large to small 
electron-hole dephasing. 

We construct the MAR solution of the Keldysh-Green's 
function equations and evaluate the dc current for arbitrary 
7 by means of the second-order recurrences similar to that for 
ballistic structures.'*'^ For the case of small electron-hole de- 
phasing (7^ 1) we show that the current is given by a general 
formula for mesoscopic connector^'* [see Eq. ( fTOl l below], i.e., 
the average of the result for the single channel'*'^ over the dis- 
tribution of transmission eigenvalues for a double-barrier nor- 
mal diffusive structure. This result coincides with the result 
of Ref. 9 and it corresponds to a broad Breit-Wigner resonance 
in the single channel. In the opposite case of large electron- 
hole dephasing, 7^1, the current can be explicitly presented 
as a sum of multiparticle tunnel currents which scale with the 
effective tunneling parameter 7^ ' . 

The structure of the paper is as follows. We formulate ba- 
sic equations and construct analytical solutions in Section HIl 
In SectionHni we present and discuss the results of numerical 
calculation of the current-voltage characteristics (IVCs). The 
multiparticle currents in the limit 7 ^ 1 are calculated ana- 
lytically in Section |IV] which also includes evaluation of the 
excess current. 



II. BASIC EQUATIONS AND SOLUTION 

We start our quantitative consideration with equation for 
the Keldysh-Green's function G(x,fi,f2) in the normal wire 

(-d <x<d). 



given in the -representation: A{E,t) = / dxe'^'^A(t + 
T/2,f — t/2). In Eqs. (|5]) we use the antisymmetric gauge 
of the superconducting phase (pR = —(pi ^ eVt, satisfying the 
Josephson relation (p = (pR — (pi = 2eVt. 

Solution of Eqs. (|2|-(|5l), being generally difficult, essen- 
tially simplifies in short junctions with opaque barriers. Av- 
eraging Eq. (O along the wire and using Eqs. (|4) and ([T]i, we 
get 



2Y[a,E,G] =/A([G^,Gk] + [G-./,Gi]) 



(6) 



where G denotes spatially averaged value of G. In the tun- 
nel limit R ^ Rn, the Keldysh-Green's functions are approx- 
imately spatially homogeneous within the normal wir o^'i^^ 

G Ri Gd ^ G-d- Thus, denoting these quantities by a single 
notation G, we arrive at the commutator equation 



[i, G] = 0, A = G+- ia.TdE, G± 



^{Gr±Gl). (7) 



A similar approach has been used in analysis of current trans- 
port in a NINIS structure. Following Refs. I27I andl28l we 
find a physically relevant solution of Eq. (|7]i satisfying the nor- 
malization condition G^ = 1 in Eq. (|2|i, 

G=^ = - r dXK{X), K{X)^{A + iX)-K (8) 

Applying Eq. (|4]i to Eq. ^ and using Eq. (H)), we symmetrize 
the quantity /(f) with respect to the left and right reservoirs, 

dX 



8eR 



TrT^[^(A),G_](f,f)- 



(9) 



The structure of the matrix current [A'(A), G_] in Eq. ^ is 
quite similar to the solution of the MAR problem for a single 



[a-E,G] ^iDd_^[Gd_^G) , G^ = l, ^= ] ■ (2) ballistic channel with the transparency D = (A^ + 1)"' given 



Here, g^'^ are the retarded/advanced Green's functions, / is 
the matrix distribution function, G^ — g^f — fg^, and the 
kernel of the energy operator E is E{ti,t2) = idt,8{ti — f2)- 
All products in Eq. ^ are time convolutions: {AB){ti,t2) = 
J dtA{t\.t)B{t,t2)- The electric current is defined as 

I{t) = {ngN/'^e)TrTK{Gd,G){t,t), Xk = o.x^ (3) 

where gAf is the conductance of the normal wire per unit 
length, and the a and T matrices operate in the Nambu and the 
Keldysh space, respectively. At the tunnel barriers x = ±.d, we 
apply the Kupriyanov-Lukichev's boundary conditions'^ 



in Ref. 15 and differs from the latter by an additional term 
—ia^XdE in the matrix K, which describes the electron-hole 
dephasing during the dwell time Xd- If this effect is negligi- 
bly small, 7^0, Eq. (|9]l rewritten in terms of the functions 
Gl.r and the transparency variable D can be transformed to 
the known formula for a short connector''* generalized to the 
nonstationary case of voltage biased SINIS junction. 



lit) = / dDTiTfc ^ , , ^ — ^ , — - 

SeRJo ''l + f({Gz„G„}-2) 



it,t), (10) 



PiD) = 



1 



gN {Gd,G)^^ = ±i2R)-' [G±d,GR,L] 



(4) 



Dp{D)dD=l. (11) 



The equilibrium Keldysh-Green's functions Grx in the right 
and left superconducting electrodes are constructed with the 
local-equilibrium Green's and distribution functions 

gR,L = ovM(e±) + ;exp(±/CTj.eyf)ayv(£'), (5a) 
./«,L=/(e±), /(£)=tanh^, e±=£±a,^, (5b) 



The fact that the MAR current in this limit is given by a con- 
volution of non-resonant single channel current with the trans- 
parency distribution p (D) for a double-barrier potential'^ (see 
also Refs. 8 and 29) is consistent with a wide resonance in the 
single channel. This result justifies the method and the result 
of Ref. 9. We note that in the static limit ^ = 0, Eq. ^ repro- 
duces the result" '^ of a direct solution of Usadel equations. 



^^■^-J(£±/0)2-A2, (5c) 



8 = 



a,E + iayAiE,(p) 
V£2-A2(£,0) ' 



(12) 
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where A(£', 0) = Acos(0 /2)[1 - iYlv{E)]-^. 

In the general case of arbitrary 7, calculation of the matrix 
^inEq. ® can be performed by expanding all quantities over 
the harmonics of eV: A{E,t) = I,„A(£',m)e-™''^'. In this 
representation, the time averaged (dc) current reads as 



/ = 



dX dE 
. \6KeR 



^ Tr^( A , £ , m) [G_ (£ , -m) , T^f ] , ( 1 3) 



and the local-equilibrium functions contain only three har- 
monics, m = 0,±1. In Eq. ( fT3] ) we rearranged the factors in 
the integrand using the fact that time averaging is equivalent to 
the trace in the time domain. After some algebra, we find the 
function G+ = ofi^ 5„,^q + iOyG\ 5\,„\ i and the commutator 
[G_,Ta:] = C7j;Gq 5m,o + /c7,,G]"m5|„,|j, where 



-N+F++N-F^ 



(14) 
(15) 
(16) 



Go = (/+ - /- A^- ) + iTy {N+ 

0+ ^^{iM + MF), Gr = im, + Mf, 

F^T, + 2/T+ , T+ = ( 1 /2) (t, + /T,) , 

A' = ReM*, M = Rev^, N = lmu^, M = Imv^ 



Here and in the following, the lower indices ± denote the en- 
ergy shift by ±ey/2. The function A^(£') istheBCS density of 
states (normalized over its value in the normal metal), which 
turns to zero at \E\ < A along with the function M{E), while 
the functions N{E) and M{E) vanish outside the energy gap. 
This leads to the following expression for the dc current. 



dEJ{E), J{E) 



dX 



-00 l6TZeR 



(17) 

j = jo + ji+j-u (18) 
;o = Tr^(£,0)(7,Go, j±i ^ ±TrK{E ,Tl)i(yyGi . (19) 

Here, J{E) is the current spectral density, whereas the quantity 
j{E,X) can be interpreted as a generalized spectral density 
depending on the auxiliary parameter X. 

According to Eq. (HJ, the matrix K{X) obeys the equation 
(A + iX)K{X) = 1; in the (E , m)-representation, it has the form 

G+[E + (m — m')eV /2,m']K{E — eVm /2,m — m) 
+ i[X - a^TdiE + eVm/2)]K{E,m) = 5,„,o 

(we omit the 'check' on top of the 4x4 matrices), where the 
sum actually contains only three nonzero terms with m' ~ 
0,±1. Introducing the quantity K,„{E) ~ K{E +meV /2,m), 
we obtain the 4x4 matrix recurrence relation 

{H,„ + iX)K,„ + h,„K„j_i + /Jm+i^Tm+i = 5,„,o, (20) 

H,,, = CT^Qm, Qm = (£■„,) — iTi/Em, 

h„, = iGyqm , qm = Gj {E,„_ 1/2), E,„ ^ E + meV. 

Solution of Eq. ( l20l i can be found by the matrix version of 
the chain fractions formalism**'^ using the ansatz 



S,„S„,-i---SiKq, ;«>0, 
P,„P,n+i-..P-iKo, m<0. 



(21) 



Recurrence relations for the "matrix chain fractions" 5„, and 
P,„ with the boundary conditions 5,„ — > at m — > 00 and P,„ 
at »n- —00 follow from Eqs. (|20] | and (l2Tl i at m 7^ 0, 



Pm = - {H„, + iX + h^Pn- 1 ) ' /!, 



m+1 • 



(22a) 
(22b) 



At m = we obtain a nonuniform equation, the solution of 
which is /To (£■ ) = (//o + /A + /lo^"- 1 + /! 1 5 1 ) " ^ . Thus the func- 
tions K{E,m) in Eq. ( fT9l ) read as 

K{E,Q)=Ko{E), K{E,l)=Si{E^)Ko{E^), (23a) 
K{E,-l)=P^i{E+)Ko{E+). (23b) 



The 4x4 recurrences in Eqs. (I22t can be reduced to the 2 x 
2 form in the Keldysh space. Indeed, assuming S,„ = —OxSm 
and P,„ = —GxPm, we arrive at the recurrences for S and P 
which are diagonal in the Nambu space, 

Sm (Qm - iXo- - qm+\S',„+\) ' Qm, (24a) 
Pm = {Qm - iX(y, -qmP'„,_iy^qm+l, (24b) 

where the prime sign denotes the change of the sign of the (7,- 
component of the matrix. Then the function Kq is also found 
to be diagonal in the Nambu space. 



(7\ -1 



(25) 



a=±l 



The functions S " and P satisfy Eqs. ( l24l i in which is re- 
placed by the scalar CT. Then, introducing the quantities Sm(A) 
and Pm{^) according to = s„,((7A) and P,^ = /7,„(c7A), we 
arrive at the 2x2 recurrences for the Keldysh matrices 



[Qm + /A ( - 1 )'" - qm+l^m+l] ' 



Pm = [2m + A ( - 1 ) "' - qmPm- 1 ] 



m+1 ■ 



(26a) 
(26b) 



In Eq. ( l25T l rewritten through the matrices 5„, and /?„,, we can 
replace aX X that does not change the result of the inte- 
gration over A in Eq. (T7\ : as the result, only the term propor- 
tional to (7. survives in Eq. 



Ko = a.K, [Qq + iX - qQP-i{X) - qisi{X) 



-1 



(27) 



By combining Eqs. ST9[ . ( l23T l, and ( l27b . and shifting the en- 
ergy in j±i by ^eV /2 which holds the result of integration 
over E in Eq. ( [TtI i unchanged, we obtain current spectral den- 
sities 

;o = 2Tr,/:(£,A)Go(£), (28a) 
ji = -2Tr, A)/:(£,A)G^(£_), (28b) 
7_i =2Tr,ii(£,A)/:(£,A)G7(£+). (28c) 

III. NUMERICAL RESULTS 

Now we proceed with the numerical analysis of the dc cur- 
rent using Eqs. (l26]l-(l28Tl and ( fTTI i. The results obtained for a 
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FIG. 1. (Color online) Current- voltage characteristics of a short dif- 
fusive SINIS junction at different values of the parameter y. 
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FIG. 2. (Color online) CuiTent vs inverse voltage in logarithmic 
scale. 



wide range of values of the parameter 7 and at zero tempera- 
ture are shown in Figs.[T]and|2] As expected, at small electron- 
hole dephasing 7= 0.2, the IVC is close to the results^ found 
from Eqs. (fTol l and (fTTT l. In this case, the IVC consists of con- 
cave portions between weakly pronounced features (steps and 
peaks) and reveals the excess current at large voltages. Since 
the edges ztA;, of the minigap « 0.8A are close to the edges 
±A of the superconducting energy gap, the existence of such 
large minigap does not distort the periodicity of the IVC fea- 
tures, the positions of which approximately coincide with the 
gap subharmonics eV = 2A/n, n = 1,2, . . .. 

As the junction transparency decreases (i.e., 7 increases), 
the excess current also decreases, approaching zero at 7sa 1, 
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FIG. 3. Current spectral density J{E) [normalized on (2ei?)^'] at 
7= 10 and different applied voltages: (a) eV = 2.5A, single-particle 
current; (b) eV = 1.5A, 2-particle current; (c) eV = 0.8A, 3-particle 
current; (d) eV = 0.6A, 4-particle current. Visible are small contri- 
butions of higher even-particle processes into the spectral density of 
the odd-particle currents. 



and then becomes negative (deficit current). Simultaneously, 
the peaks in the IVC almost vanish, and the subharmonic gap 
structure at small eV and 7 1 — 3 becomes somewhat chaotic 
(see Fig. |2]i. This is due to the interplay of the contributions 
of MAR trajectories touching the superconducting gap edges 
and the edges of the minigap Ag « (0.25 — 0.5)A. Then, at 
large 7=10 — 30, the IVC features become regular again and 
their positions exactly correspond to the gap subharmonics. 
In this case, the minigap is small, A^ ss 0.03 — 0.1, and there- 
fore affects the MAR trajectories touching the superconduct- 
ing gap edges with even number of steps only. The enhanced 
density of states in the vicinity of the minigap increases the 
transmissivity of these MAR chains; this leads to anomalous 
enhancement of the magnitude of the dc current just above the 
even gap subharmonics n = 2k. As k increases, this resonance 
effect becomes more pronounced and leads to the appearance 
of the IVC portions with negative differential resistance, as 
seen in Fig. |2] The current spectral density J{E) shown in 
Fig.|3]at 2A/n < eV < 2A/(n - 1), n = 1 -4, has the form of 
n main equal peaks which acquire a resonant shape for n = 2k 
[Fig EIb) and|3ld)], in accordance with the above-mentioned 
anomalous transmissivity of even MAR chains. Small foot- 
prints of these resonances are visible in the spectral density of 
the odd-particle currents with n = 2k—\ [Fig. [3a) and|3lc)]. 

We note that, at 7 > 1 and eV < A, the averaged IVC is well 
approximated by the dependence 



/(y) 



O.ISA^ 
"eRAT 



-2A/eV 



(29) 
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which is similar to the resuh for the ballistic SIS structure^" 
with 7^' standing for the transparency D. Thus, at 7 > 1, 
the quantity 7^ ^ plays the role of an effective parameter for 
multiparticle tunneling processes which determines the value 
of the n-particle current of the order of R^^y^^" . This con- 
clusion is confirmed by asymptotic analysis of multiparticle 
currents (see below). 



MaMp with /3 7^ a), which turn to zero at eV °°, 



X 1 + 



Zk 

3/2 



2M_ 



1/2^1/2.-3/2- 



-1/2 



+ ^-1/2.-3/2 1 1 - 



1/2 
1/2 



-1/2,-1/2 



(36) 



-1/2 



-4- 



Lal}=NaNii{fa-fp): a > p . 



-1/2,-1/2 



(37) 



IV. ANALYTICAL RESULTS 

In this section, we present a detailed analytical discussion 
of the dc current. To this end, we consider Eqs. (l26t as func- 
tional equations for the functions S{E,X) — siG'l{E+) and 
P(£,A)=/?_iG+(£_), 

{S,P){E,X)= [K^\e ±eV,-X) +Yl^'''{E ±eV,-X)Y] (30) 

(31) 



K^\E,X)=G^-iTdE + iX, 
n^-^(£,A) = -Gi{E±){S,P){E,X)Gi{E±). 



(32) 



In terms of these functions, the contributions Eq. (I28bt and 
(12 Sel l of nonzero harmonics to the dc current read as 

;±i = T2Tvr{P,S)GiiE^)K-'G^{E^), (33) 

(34) 



K ^k'-^ +n^. 



Analysis of Eq. (l30l l shows that the matrices /T", P, and 
S possess certain symmetry properties with respect to the 
transformation E -E, X -A: ^ -5''+, P- ^ S\ 
Zp o Zs, Zk ^ Zk, where Zk = detK, Zp = detP, Z5 = det5, 
and the upper indices 1, z, and + denote 1-, t^-, and T+ -matrix 
components, respectively. These relations allow us to exclude 
the antisymmetric terms that vanish under integration over E 
and X in Eq. ( [TtI i and to write down the spectral densities 
Eqs. (l28T l in a compact form 



Here and in the following, the lower indices denote the en- 
ergy shift, e.g., Na = N{E + aeV). At zero temperature 
(/ = sgnii), the factor L^p is nonzero within the energy re- 
gion 



A-aeV <E < -A-peV. 



(38) 



Existence of this energy interval, in which Lap = 2NaNp, im- 
poses the following condition on the applied voltage. 



{a-p)eV >2A; 



(39) 



otherwise, the function L^p is identically zero. Since the con- 
vergence of integration over E of all terms in curly brackets in 
Eq. d36T l is ensured by the functions Ma or Ma with Ea ^ A, 
the functions Np with different energies (j3 a) can be ap- 
proximated by their values in a normal metal Np — 1. By 
using these simplifications, we arrive at the following asymp- 
totics of the dc current. 



2A f°° reV/2A ri 

/=-^y dX dxji{x,X) + J^dxj2{x,X) 

x[{X+y)^+z'^-l] _. 2yi 



(40) 



|A2-l + (z-/3;)2|2' 

= 2yx\J — 1, z = X ^ 



n = 



|A2 + l-(zi-,3;i)2|2' 



= VT^, Zi =x(l+27yi). 



By integrating over X in Eq. ( |40] | and separating out the con- 
stant term, we obtain a general expression for the excess cur- 
rent. 



./o = -—\2K\N+U -N-f-)-K+{N+ -N-)] , (35a) 



j,+U = -^M^{P^K+-P+K^) 



A. Excess current at eV 2> A 



(35b) 



r-(7) = 



h [A 



eR 

y/2x(T + T 



/ dx j\ [x) + dx j2{x) — 1 
J\ ' Jo 



(41) 



tVt-a 



1, jlix) 



2^2(1- 



T^=A\+Ay\zl Ar=l+y\~zl 



We start with evaluation of the excess current I"^'^ at large 
applied voltage and for arbitrary 7. This quantity is con- 
tributed by both the single-particle current and the two-par- 
ticle Andreev current. Formally, Z*"" is the voltage-indepen- 
dent term in asymptotic expression for the dc current / = 
V /2R + 1'''"' + 0{A/eV) at eV > A. In coiTesponding expan- 
sion of full current spectral density Eq. (|35t . we truncate the 
recurrences for P and S, omitting in j{E,X) the combinations 
of the functions Ma or Ma with different energies (such as 



In the limit of small dephasing 7 — > 0, when the integral 
overx in Eq. (|40] | can be explicitly calculated, the substitution 



D= (A 

rexc 



2 -I- 1 ) ' leads to the formula^ 

-1 



^J^dDpiD)IlfsiD) ^0.53^, 7«1, (42) 



IliliD) 



D^A 



£)2 



2(l+3^)V^"'l-V^ 



In 



, 3?= 1 - Z), 
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which expresses the excess current through its value for a 
single ballistic channel^' averaged over the transparency dis- 
tribution Eq. ( fTTT l. in accordance with Eq. (fTOl i. In the opposite 
case 7^1, becomes negative (deficit current), 



3eR' 



7> 1. 



(43) 



determined by the last term in Eq. (l46T l is similar to that for the 
function L^p [see Eqs. (l37T i and (I38]l1. therefore this term turns 
to zero at eV < 2A, according to Eq. ( |39] |. Thus, in this voltage 
region, the action of matrix envelopes G+o on the matrix 
is reduced to multiplication on the scalar factor — M+. Similar 
considerations applied to each term of the expansion Eq. (145 al) 
lead to the following simplified series for the functions fl. 



B. Multiparticle currents at large 7 

In the limit of large dephasing 7^ 1, it is possible to ex- 
press analytically the full current as a sum of contributions 
of n-particle tunneling processes.^ Here, we proceed with the 
asymptotic analysis of these partial contributions. First, we 
separate out the unity and the traceless components of the 
Keldysh matrices, e.g., P = +P, P= T-^P' + 



n^'^ = -G±o(5,P)G±o, n^''^^=M±i/2(5',P^), (44a) 

,2 1/^2 .,2\ 



Ma 



[g+(£„)]' = -(m;-m^,, 



(44b) 



where G±„, — G\{E± ± meV). By introducing the notations 

Z±„, EE det{^!^„,+7i:0;„ + n^'^[£±,„,(-l)"'A]} , m > 0, 

we rewrite Eqs. (|30] l in an expanded form, explicitly perform- 
ing the recurrences for the functions P and S. The result can 
be presented in the form of the series for the functions 11, 



= til n^' = L </ ^A^a-l/2z„^ (45a) 

m=l m=l o:=l 

m 

= G+oGi...G„,_iX,„G,„_i...GiG+o ]~[ (45b) 

a=l 

The series for 11^ differ from Eqs. (l45i by opposite signs of all 
lower indices. With this remark, the series for the functions P 
and S can be obtained from Eqs. (l45T l and (|44] |. These series 
can be interpreted as asymptotic expansions over 7^ ^ , due to 
the presence of large parameter 7^ 1 in Z^. Physically, these 
expansions reflect the nature of the net current as a sum of 
n-particle tunnel currents;^ each of them exists at eV > 2A/n 
and scales as 7^^" with respect to the single-particle current. 
The latter fact allows us to consider the n-particle current only 
within its actual voltage region 2A/n < eV < 2A/(n — 1); at 
larger voltages, the (n — l)-particle current dominates. Es- 
timation shows that mth terms in Eqs. (l45l l contribute to the 
{m+ l)-particle current; thus, it is enough to consider them 
only at eV < 2A/m, which greatly simplifies the structure of 
the series. 

Indeed, consider the term flj proportional to the product 



G+o^:l'G+o = -M+K'; - -NaM'p (fa -fp)T+ 



(46) 



where a — 3/2, P — 1/2, and we used the identities = I 
and FaFpFa = + {^fa — 2fp ) T+ . The allowed energy region 



tl'^Y. ei2A~meV){~iy"Kl!,YlMa,-l/2Za 



(47) 



where we introduced the Heaviside step function to spec- 
ify explicitly the relevant voltage regions. The series for the 
functions fl^ differ from Eq. (|47] i by opposite signs of lower 
indices. 

Now we proceed with asymptotic evaluation of the dc cur- 
rent. First we consider the contribution 70 in Eq. (|35a| l to the 
net current spectral density which, according to Eq. (|34i , can 
be presented as a sum of three terms, 



jo- .10+ Jo + Jo, 



(48) 



;,f = - — [2K^ {N+U -N-f-)- K+ {N+ -N^)]. (49) 

The first term is equal to — (4/Za:)-Li/2,-i/2 ^nd represents 
the spectral density of the single-particle current. Accord- 
ing to Eqs. ( |38] ) and ( [39l ), it exists within the energy interval 
\E\ < —A + eV/2 and vanishes at eV < 2A. Thus, at subgap 

PS 

voltages, we have to involve the terms Jq , which differ from 
Eq. ( |49] l by replacements Kq — ?> n^"^. Considering, e.g., the 
spectral density /q and taking into account Eq. (|47] |. we found 
that the mth term in the expansion of Jq is proportional to 



9i2A~meV) £ /3L,„+„,j 

a,j3=±l/2 



(50) 



As follows from Eq. ( [39] l, only the term with a = — j3 = 1 /2, 
proportional to 9 [{m + 1 )eV — 2A] , survives in Eq. ( fSOl l. Thus, 
we obtain the following series for j^. 



Jo 



'a-l/2 

.y /v"i\- /-^m-1/2,-1/2 X \7~\ 



m— 1 

n 



, «a = ^, (51) 



1, 2A/m<ey <2A/(m-l), 
otherwise. 



By applying similar considerations to and i\ + we ar- 
rive at the formula for full generalized current spectral density 
[at n = 1, the product in Eq. (|53] l is unity], 



y(£,A) = £;t„(y)/"), 



4 " 



IZ^I 



m—X 

n 



«a-l/2 



(52) 



(53) 



According to Eq. (|53] l, the n-particle spectral density /"^ 
consists of n equal contributions of MAR chains with n steps. 
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Each chain starts at the energy 1/2 < —A and finishes 

at > A> thus transferring the quasiparticles to the ex- 

tended states above the energy gap, which results in formation 
of the dissipative current. The intermediate energies /ia-1/2 
inside the gap correspond to the points of Andreev reflections. 
These contributions are nonzero within the energy intervals 
of width neV — 2A, which are distributed equidistantly (with 
spacing eV) along the energy axis and symmetrically with re- 
spect to the zero energy, in conformity with the numerical re- 
sults shown in Fig. [3] This enables us to write down the full 
dc current as a sum of n-particle tunnel currents /("\ where 
only one term in multiplied by n is taken into account. 



(n) 



(54) 



«=i 



-00 In jK-eV/l 



eR \Zk\ 



a=\ 



(55) 



At « = 1, the product in Eq. (I55l l is assumed to be unity. 

To complete our consideration, we present final expressions 
of n-particle currents for n — l,2, and 3 obtained from Eq. (ISSl l 
by the integration over X. A nontrivial point in this procedure 
is a proper choice of approximation for the determinants Z. In 
the single-particle current, it is enough to take Zk in the main 
approximation as detK° = -[{X- TdE f + {1 /4) (A^+ + A^- )^] , 
neglecting contributions of n^"^ to the function K in Eq. ( l34b . 
As is obvious from this expression, the parameter drops out 
from and we obtain a simple formula 



dE 



-f eR{N_ 



(56) 



From the standpoint of the circuit theory for incoherent SINIS 



structures,'^ this result can be interpreted as the Ohm's law 



for two tunnel resistors R± = RN^^ connected in series. 

In calculation of the two-particle current, the main ap- 
proximation is applicable to the determinant Z i = — 
whereas in Z/f = — \Xk'\^ one should hold the term 11^, 



1 

x_i = a_i + -A?_3/2, 
1- 



XK^iXK + ]rN+ + ^, 



Xk^X-t^ciE^-N-, X-x=-XK-2XdE-+N-. 

By taking Xk as a new integration variable, we see that its 
characteristic value is of the order of unity, which enables us 
to approximate A_i as — 2T(/£'_. After integration over Xk and 
symmetrization of the allowed energy interval, we obtain 



,(2) 



dE 



A-eV eR Ni [(4Trf£)2 +Nli] +4aN^i 



(57) 



The current spectral density in Eq. (157) has a resonant form, 
with a sharp peak at zero energy (the resonant nature of the 
even-particle currents has been already noted in Sec. HITt . If 
the applied voltage is not very close to the threshold A/e of 



the two-particle current, the integral in Eq. (l57i can be calcu- 
lated in the resonant approximation by assuming £ = in all 
spectral functions and spreading the limits to ±00, 



^(2) 



TtA N{eV) 



(58) 



Similar considerations lead to the following expression for the 
three-particle current. 



/(3) 



Ay^eR 



-A+leV dE NiO+a^N 3 
2 



A-^eV a^N ^Ei+a+NiEi 



(59) 



Expressions Eqs. (|56]|-(|59]| well reproduce the results of full 
numerical calculations for large 7 shown in Figs. [T| and |2l 



V. SUMMARY 

We have calculated, both numerically and analytically, the 
current-voltage characteristics (IVCs) of a diffusive SINIS 
junction, where S are local-equilibrium superconducting re- 
servoirs, I denotes tunnel barriers, and N is a short diffusive 
normal wire, the length of which is much smaller than the 
coherence length and the resistance R^ is much smaller than 
the resistance R of the tunnel barriers. The regime of co- 
herent MAR transport in such structure is governed by the 
parameter 7 — T^/A, which represents a characteristic phase 
shift between the wavefunctions of the electron and the retro- 
reflected hole accumulated during the quasiparticle dwell time 
Xci ^ E^^ (R/Rn). We demonstrated that the Keldysh-Green's 
function equations for this problem can be efficiently solved 
in the whole range of electron-hole dephasing < 7 < 00. This 
is achieved by reducing the solution of full 4x4 matrix two- 
time Keldysh-Green's function equations'^ to the solution of 
the 2x2 matrix recurrence relations of the second order, sim- 
ilar to the recurrences in analogous ballistic problems. In 
the limit of small dephasing 7^0, our solution reduces to a 
known formula for mesoscopic connector, i.e., averaging of 
the result for the single channel junction over the distribution 
of transparencies for the corresponding double-barrier normal 
diffusive structure.^ 

In the opposite case of large electron-hole dephasing 7^1, 
the subharmonic gap structure in the IVC scales with 7^'; this 
means that the n-particle tunnel currents scale as 7'^" with 
respect to the single-particle current, and 7^' plays the role 
of an effective tunneling parameter Due to the presence of 
resonant MAR chains touching the edges of small minigap 
Ag « A/ (7+ 1 ), the even subharmonics are enhanced, and cor- 
responding portions of the IVC show negative differential re- 
sistance. We presented analytical results for the excess current 
at arbitrary 7 and for multiparticle currents at 7 1 . 

For experimental observation of the phenomena discussed 
in this paper, the most stringent constraint concerns the Jo- 
sephson regime. This implies small values of the dwell time 
compared to the inelastic relaxation time. If this require- 
ment is not fulfilled the central metallic island acts as a reser- 
voir, and the structure splits in two NIS junctions connected 
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in series. This is the case of SINIS junctions extensively 
used in microcoolers?? and SET turnstiles. For the Joseph- 
son effect to occur in metallic SINIS junctions with conven- 
tional oxide tunnel barriers, a sandwich-type junctions must 
be employed having extremely thin normal metallic layer not 
exceeding 10 nm. Such junctions have been developed us- 
ing Nb/AlO.r/Al/AlO^/Nb technology, and they demonstrated 



rather large values of 7 ~ 10 and a pronounced deficit 
current. ^"""^^ This is precisely the limit of large electron-hole 
dephasing studied in this paper In order to investigate a 
crossover to the regime of small dephasing at 7 ^ 1 one 
needs to use junctions with more transparent NS interfaces, 
such as junctions based on diffusive InAs 2D electron gas or 
graphene, or corresponding nanowires and nanotubes. 
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